/*
 * UserDefinedFunc.cpp
 *
 *  Created on: Nov 11, 2011
 *      Author: yi
 */
#include "mbuild/UserDefinedFunc.h"


ChargeDensity::ChargeDensity()
{
	this->nlayer=0;
	this->density=NULL;
	this->offset=NULL;
	this->std=NULL;
	scaleFactor=1;
	margin=1e-10;
	IsAdditive=true;
}



ChargeDensity::ChargeDensity(int nlayer)
{
	this->nlayer=nlayer;
	this->density=new double[nlayer];
	this->offset=new double[nlayer];
	this->std=new double[nlayer];
	scaleFactor=1;
	margin=1e-10;
}



 ChargeDensity::~ChargeDensity(void)
{
	 delete [] density;
	 delete [] offset;
	 delete [] std;
}
 void ChargeDensity::set_nlayer(int n)
 {
 	nlayer=n;
 }


void ChargeDensity::set_layer(int layer, double offset, double width, double sigma)
{
	/*
	 * set up layer i with offset,width,density.
	 * layer: 0->nlayer
	 */
	set_width(layer,width);
	set_sigma(layer,sigma);
	set_offset(layer,offset);
}



void ChargeDensity::set_width(int layer, double width)
{
	if(width <=0) throw ParameterError();
	(this->std)[layer]=width;
}



void ChargeDensity::set_sigma(int layer, double sigma)
{
	(this->density)[layer]=sigma;
}



void ChargeDensity::set_offset(int layer, double offset)
{
	(this->offset)[layer]=offset;
}



void ChargeDensity::skipMargin(double margin)
{
	if(margin >1) printf("Warning: Charge Density Margin too large!\n");
	this->margin=margin;
}



bool ChargeDensity::IfEqualSolute(double v)
{
	if(fabs(v)<margin) return false;
	return true;
}



void ChargeDensity::set_density(double f)
{
	scaleFactor=f;
}

double ChargeDensity::get_density()
{
	return scaleFactor;
}

double ChargeDensity::f(double d)
{
	double v=0.0,d0;
	for(register int i=0;i<nlayer;i++) {
		d0=(d-offset[i])/std[i];
		v+=scaleFactor*density[i]*exp(-d0*d0/2)/(2.50662827*std[i]);
	}
	return v;
}


DielectricConstant::DielectricConstant()
{
	in=2;
	out=80;
	on=10;
	width=3;
	mode=SOLUTE_EQL;
	IsAdditive=false;
}

DielectricConstant::~DielectricConstant()
{
}



void DielectricConstant::set_headgroup_width(double w)
{
	if(w<0) throw ParameterError();
	width=w;
}



void DielectricConstant::set_headgroup_diel(double hg)
{
	if(hg<0) throw ParameterError();
	on=hg;
}



void DielectricConstant::set_water_diel(double water)
{
	if(water<0) throw ParameterError();
	out=water;
}



void DielectricConstant::set_memb_diel(double memb)
{
	if(memb<0) throw ParameterError();
	in=memb;
}


void DielectricConstant::set_solute_diel(double solute)
{
	if(solute<0) throw ParameterError();
	solute_value=solute;
}


void DielectricConstant::set_diels(double water, double hg, double memb,double solute)
{
	set_water_diel(water);
	set_memb_diel(memb);
	set_headgroup_diel(hg);
	if(solute == 0) solute=memb;
	set_solute_diel(solute);
}



Kappa::Kappa()
{
	solute_value=0;
	IsAdditive=false;
	in=0;
	out=1;
	on=0;
	width=3;
	mode=SOLUTE_EQL;
}



Kappa::~Kappa(void)
{
}


Kappa::Kappa(double width)
{
	solute_value=0;
	IsAdditive=false;
	in=0;
	out=1;
	on=0;
	mode=SOLUTE_EQL;
	set_headgroup_width(width);
}

void Kappa::set_headgroup_width(double w)
{
	if(w<0) throw ParameterError();
	width=w;
}

IMM1Transition::IMM1Transition()
{
	in=0;
	out=1;
	nsmth=10;
	thick=27;
	mode=DONOTHING;
}



IMM1Transition::IMM1Transition(int nsmth, double thick)
{
	in=0;
	out=1;
	mode=DONOTHING;
}



IMM1Transition::~IMM1Transition()
{
}



void IMM1Transition::set_memb_thickness(double thick)
{
	if(thick <=0) throw ParameterError();
	this->thick=thick;
}



void IMM1Transition::set_nsmth(int nsmth)
{
	if(nsmth <=0) throw ParameterError();
	this->nsmth=nsmth;
}



void IMM1Transition::set_transition_width(double w)
{
}

GaussianGeoConstraint::GaussianGeoConstraint()
{
	mode=DONOTHING;
}



GaussianGeoConstraint::~GaussianGeoConstraint()
{
}



GaussianGeoConstraint::GaussianGeoConstraint(double forceconstant, double offset, double width)
{
	set_forceconstant(forceconstant);
	set_offset(offset);
	set_width(width);
}



void GaussianGeoConstraint::set_forceconstant(double f)
{
	this->forceconstant=f;
}



void GaussianGeoConstraint::set_offset(double o)
{
	this->offset=o;
}



void GaussianGeoConstraint::set_width(double w)
{
	if(w<=0) throw ParameterError();
	this->std=w;
}



double GaussianGeoConstraint::f(double d)
{
	d-=offset;
	d/=std;
	return (1-exp(-d*d))*forceconstant;
}



SwitchGeoConstraint::SwitchGeoConstraint()
{
	mode=DONOTHING;
}



SwitchGeoConstraint::~SwitchGeoConstraint()
{
}



SwitchGeoConstraint::SwitchGeoConstraint(double forceconstant, double cutoff)
{
	set_forceconstant(forceconstant);
	set_cutoff(cutoff);
}



void SwitchGeoConstraint::set_forceconstant(double forceconstant)
{
	in=forceconstant;
	out=forceconstant;
}



void SwitchGeoConstraint::set_cutoff(double c)
{
	width=c;
}



HarmonicGeoConstraint::HarmonicGeoConstraint()
{
	mode=DONOTHING;
	stype=BOTH;
}



HarmonicGeoConstraint::~HarmonicGeoConstraint()
{
}

HarmonicGeoConstraint::HarmonicGeoConstraint(double f,double o,double cutoff)
{
	mode=DONOTHING;
	stype=BOTH;
}


void HarmonicGeoConstraint::set_forceconstant(double f)
{
	forceconstant=f;
}

void HarmonicGeoConstraint::setmode(SYMETRICTYPE s) 
{
    stype=s;
}

void HarmonicGeoConstraint::set_offset(double o)
{
	offset=o;
}



void HarmonicGeoConstraint::set_cutoff(double cut)
{
	cutoff=cut;
}



double HarmonicGeoConstraint::f(double d)
{
	d-=offset;
	switch(stype) {
	case BOTH :  d=fabs(d)-cutoff; break;
	case INSIDE : d= d>0? d-cutoff:0; break;
	case OUTSIDE : d= d< 0?d-cutoff:0; break;
	}

	if(d>100) return (forceconstant*1e4);
	if (d<=0) return (0.0);
	else return (forceconstant*d*d);
}

QuadraticChargeDistribution::QuadraticChargeDistribution()
{
	a=0;
	b=0;
}



QuadraticChargeDistribution::~QuadraticChargeDistribution()
{
}



void QuadraticChargeDistribution::modify(Dfunc *df,double x,double y,double z)
{
	//f->apply_filter(this);
	double scalefactor=f(x,y,z);
	df->set_density(scalefactor);
}

void QuadraticChargeDistribution::setup(double a,double b,double R0,double K0,double T) {
	this->a=a;
	this->b=b;
	this->R0=R0;
	this->K0=K0;
	this->T=T;
}

double QuadraticChargeDistribution::f(double x, double y, double z)
{
	double r=x*x+y*y;
	if(r<(R0+K0)*(R0+K0))	return a+b*z*z;
	else return 1;
}

void DielectricConstant::set_density(double d)
{
	return;
}



double DielectricConstant::get_density()
{
	return on;
}



void Kappa::set_density(double d)
{
}



double Kappa::get_density()
{
	return 0.0;
}



void IMM1Transition::set_density(double d)
{
}



double IMM1Transition::get_density()
{
	return 0.5;
}



void GaussianGeoConstraint::set_density(double d)
{
	forceconstant=d;
}



double GaussianGeoConstraint::get_density()
{
	return forceconstant;
}



void SwitchGeoConstraint::set_density(double d)
{
   in=d;
   out=d;
}



double SwitchGeoConstraint::get_density()
{
	return in;
}



void HarmonicGeoConstraint::set_density(double d)
{
	forceconstant=d;
}



double HarmonicGeoConstraint::get_density()
{
	return forceconstant;
}
















